Polar meron-antimeron networks in strained and twisted bilayers

Out-of-plane polar domain structures have recently been discovered in strained and twisted bilayers of inversion symmetry broken systems such as hexagonal boron nitride. Here we show that this symmetry breaking also gives rise to an in-plane component of polarization, and the form of the total polarization is determined purely from symmetry considerations. The in-plane component of the polarization makes the polar domains in strained and twisted bilayers topologically non-trivial, forming a network of merons and antimerons (half-skyrmions and half-antiskyrmions). For twisted systems, the merons are of Bloch type whereas for strained systems they are of Néel type. We propose that the polar domains in strained or twisted bilayers may serve as a platform for exploring topological physics in layered materials and discuss how control over topological phases and phase transitions may be achieved in such systems.

Out-of-plane polar domain structures have recently been discovered in strained and twisted bilayers of inversion symmetry broken systems such as hexagonal boron nitride. Here we show that this symmetry breaking also gives rise to an in-plane component of polarization, and the form of the total polarization is determined purely from symmetry considerations. The in-plane component of the polarization makes the polar domains in strained and twisted bilayers topologically non-trivial, forming a network of merons and antimerons (half-skyrmions and half-antiskyrmions). For twisted systems, the merons are of Bloch type whereas for strained systems they are of Néel type. We propose that the polar domains in strained or twisted bilayers may serve as a platform for exploring topological physics in layered materials and discuss how control over topological phases and phase transitions may be achieved in such systems.
Recently it has been realized that ferroelectricity can occur in layered systems comprised of stacks of two-dimensional (2D) materials such as hexagonal boron nitride (hBN) (see Fig. 1a), provided the stack of layers does not have inversion symmetry 1 . In an aligned stack of hBN (3R stacking), which has four non-orthogonal mirror planes and is therefore non-centrosymmetric but still non-polar, sliding one layer over the other breaks the mirror symmetry about the plane which is parallel to and half-way between the layers, resulting in an interlayer transfer of electronic charge and an out-of-plane polarization 1-3 (see Fig. 1b). For anti-aligned hBN (2H stacking), there is an inversion center for every stacking, and the system is nonpolar. Applying an electric field to aligned hBN, the polarization can be inverted via a relative sliding between the layers (van der Waals sliding) 4,5 in order to align the polarization with the field (see Fig. 1c). This mechanism is highly unconventional when compared to the ferroelectricity observed in ABO 3 oxide perovskites, in particular, because the polarization generated is perpendicular to the atomic motion.
In a twisted bilayer, two layers are twisted with respect to one another, forming a supercell known as a moiré superlattice (see Fig. 1d-f). A moiré superlattice can also be generated by introducing a small relative strain or lattice mismatch between the layers. Twisting has been shown to result in novel phenomena such as superconducting 6 and insulating 7 behavior in bilayer 'magic angle graphene', and recently ferroelectricity in hBN 5,8,9 . A small lattice mismatch moiré superlattice formed between non-Bravais lattice monolayers has local regions with different stacking configurations, which may locally break mirror symmetry. This symmetry breaking in conjugation with the absence of an inversion center in the monolayer leads to local out-of-plane polarization with stacking-dependent direction 2,3 . Thus, the stacking domains in strained and twisted bilayers can be identified as out-of-plane 'moiré polar domains' (MPDs). The experimentally observed ferroelectricity has been attributed to the motion of the domain walls separating the MPDs in response to an applied out-ofplane electric field. As a result, the MPDs with polarization (anti-) aligned to the field (shrink) grow in size.
Something which to our knowledge has not been considered is the possibility of an in-plane polarization, both in moiré superlattices and commensurate layered systems. Because layered systems are periodic in the in-plane directions, the in-plane polarization is a latticevalued quantity, and only changes in the in-plane polarization are well defined, modulo a quantum of polarization 10 . 2D honeycomb compounds with an AB sublattice structure have a triangular in-plane polarization lattice 11 , and it is natural to expect that changing the stacking configuration in a bilayer may result in a continuous change in the in-plane polarization.
We show with first-principles calculations of bilayer hBN that an in-plane polarization is indeed generated in the layered system when one layer slides over the other. As a consequence, the MPDs do not just point in the out-of-plane direction, but also have an intricate in-plane component, such that the polarization vector has topologically nontrivial winding, and the MPDs form a network of merons and antimerons (winding numbers ± 1 2 ). This indicates that the polar properties of layered systems can have rich topological structures. Topological polar structures such as skyrmions 12,13 and merons 14 have been observed in ferroelectric materials such as oxide perovskites, and have been shown to result in novel physics of interest for future applications in nanotechnology, such as negative capacitance 15 and high-density information processing 13 . So far, band topology in the moiré systems has appeared in the electronic structure of magic angle graphene 16,17 , Chern bands in twisted topological insulators 18 , and topological superconductivity in twisted cuprates 19 . The polar meron-antimeron network suggests that moiré materials also exhibit real space topology, echoing recent similar discoveries of topologically nontrivial strain fields in twisted bilayers 20 and magnetic textures in moiré patterned topological insulators 21 .

Results
We first calculate the out-of-plane and in-plane polarization of 3Rstacked bilayer hBN, using the SIESTA 22 and ABINIT 23 first-principles codes (see the "Methods" section). The polarization is calculated in commensurate 3R-stacked bilayer hBN as a function of relative displacement s between the layers, i.e. in configuration space 24 , which can be used to estimate the polarization in real space for arbitrary strains and twist angles, provided the mismatch is small enough that the local stackings are well approximated by a commensurate bilayer plus a relative translation (see Supplementary Material, Section I). A changing in-plane polarization was found, of the same order of magnitude as the out-of-plane polarization (see Fig. 2a and b).
The shape of the polarization field as a function of relative stacking is determined purely from symmetry considerations, although the magnitude is material specific. The aligned AA stacking with space group P 6m2 (#187) has three out-of-plane mirror planes running through C 3 rotations of thex +ŷ unit cell diagonal, wherê , and an in-plane mirror plane halfway between the layers. The bilayer is therefore non-polar for this stacking, but because the mirror planes are not orthogonal, it is not centrosymmetric. Sliding one layer over the other by 1 3 or 2 3 along thex +ŷ unit cell diagonal (or one of its C 3 rotations), the energetically favorable AB and BA stackings are realized, both with the polar space group P3m1 (#156). The three mirror symmetries through C 3 rotations of the unit cell diagonal are preserved, but the in-plane mirror symmetry is broken, allowing for polarization only in the out-of-plane direction. Halfway between the AB and BA stacking configurations, at the saddle point (SP, x = 1 2 ), the Abm2 (#39) space group is realized, which only has an out-of-plane mirror symmetry through thex +ŷ unit cell diagonal, assuming the relative translation of the layers is along this diagonal. Additionally, the system is left invariant after mirroring about the plane half-way between the layers plus a non-symmorphic translation of 1 2x +ŷ À Á , preventing any out-of-plane polarization. Thus, only an inplane polarization along thex +ŷ unit cell diagonal is allowed. For any other translation along thex +ŷ unit cell diagonal or one of its C 3 rotations, the bilayer has the Cm (#8) space group, with only the mirror plane running through that diagonal. The polarization is then confined to that mirror plane but can have both in-plane and out-of-plane components. Finally, for a translation not along the unit cell diagonal, the P1 (#1) space group is realized, and the polarization can point in any direction.
Using the in-plane mirror symmetry of the AA stacking configuration, we can further deduce that the out-of-plane polarization P ⊥ must be an odd function of in-plane translations (see Supplementary Material, Section I). Requiring also that P ⊥ transforms as a scalar field with respect to C 3 rotations about the out-of-plane axes through AA, b Illustration of the interlayer charge transfer Δq(x) and resulting polarization arising from relative sliding x in bilayer hBN, where q B and q N are the charges of the B and N atoms, respectively, and d 0 (x) is the equilibrium layer separation. c Out-ofplane polarization P ⊥ and d stacking energy V stack in bilayer hBN as a function of relative stacking x along the unit cell diagonal, from first-principles calculations in ref. 3 . The stacking configurations at the extrema, AA, AB, SP, and BA are labeled and sketched above the plots, and their stacking energies V AA , V SP , and V AB = V BA are marked on the vertical axis of (d). Changing between the energetically stable AB and BA stacking configurations via vdW sliding inverts the polarization between the minimum and maximum values P min and P max . e-g Sketch of red and blue hexagonal bilayers with relative twist angles of, θ = 0°, 5°, 10°, respectively. AB, and BA, it can be shown to be of the form P ? ðx,yÞ = P odd 1 sinð2πxÞ + sinð2πyÞ À sinð2πðx + yÞÞ ½ ð1Þ in fractional coordinates x in configuration space, where (x, y) ≡ x are fractions of the lattice vectorsx andŷ; the fractional and Cartesian coordinates in configuration space are related via s = gx, where We can also show that the in-plane polarization P ∥ , must be even with respect to in-plane translations. Additionally, it should transform like a vector with respect to C 3 rotations about the out-of-plane axes, and therefore must be of the form ΔP k ðx,yÞ = P even 1 cosð2πxÞ À cosð2πðx + yÞÞ Eqs. (1) and (2) were fit to the 1D data along thex +ŷ unit cell diagonal in Fig. 2a and b. The 1D data were sufficient to obtain a good fit for the 2D functions, which was verified by fitting to a larger set of displacements parametrizing the entire unit cell in configuration space. The out-of-plane and in-plane polarization in Cartesian coordinates, are obtained by the transformations P ⊥ (s) = P ⊥ (g −1 x) and ΔP k ðsÞ = g À1 T ΔP k ðg À1 xÞ and are shown in Fig. 2c and d, respectively.
The shape of the out-of-plane polarization in bilayer hBN as a function of relative stacking is well known 2,3 : P ⊥ (s) forms a triangular domain structure, with each domain having three neighboring domains of opposite polarization. The in-plane polarization has a remarkable structure: ΔP ∥ (s) flows into and out of the centers of the AB and BA domains, at which it is zero. The magnitude of ΔP ∥ (s) is maximal along the lines joining the centers of the AB and BA domains. A six-pointed star forms around the AA stacking configuration, at which both inplane and out-of-plane components of the polarization are zero, modulo a quantum of polarization. Also, we note from symmetry that ΔP ∥ (s) ∝ ∇ s P ⊥ (s).
Mapping from configuration space to real space, we see that the total polarization has a different form, depending on whether the moiré superlattice is induced via relative straining or twisting. For a relative strain η between the layers, the mapping between configuration space and real space is s = ηr. In this case, configuration space and real space are related via a simple scaling by η. In Fig. 3a and b we show the out-of-plane and in-plane polarization in a strained bilayer, which is identical to the polarization in configuration space. For a relative twist θ, the mapping between configuration space and real space is For scalar fields, the quantities are related via scaling by θ and a reorientation of the cell vectors. For vector fields, the general form is different in both spaces. In Fig. 3d and e we show the out-of-plane and inplane polarization in a twisted bilayer. We can see that the out-of-plane polarization has the same form, but the in-plane polarization curls in a clockwise manner around the centers of the AB/BA domains. The components of polarization in Fig. 2a and b were reproduced by integrating the dynamical charges, Z * κ,αβ = V ∂P α ∂s κ,β 25 (see Supplemen-tary Material, Section V). This further allows the decomposition of the polarization into the contributions from the displacements of different atoms and in different directions. The polarization is mostly generated by in-plane sliding, with negligible contributions from the out-of-plane displacements, i.e. the rippling of the interlayer separation as one layer slides over the other. As a result, the antisymmetric part of Z * (s) makes a significant contribution to the polarization, which is highly unusual for a ferroelectric material. The in-plane component suggests that the polarization field in twisted bilayers does not just point in the out-ofplane direction, but exhibits intricate winding which is topologically nontrivial. Topology has played a manifest role in 2D materials, ranging from band theory to skymrions in magnetic systems. Such skymrions arise due to a mapping from a periodic unit cell to a classifying space that is topologically a sphere, as quantified via a homotopy π 2 ðS 2 Þ = Z winding.
In Fig. 3c and f, we show the local topological charge of strained and twisted bilayer hBN, respectively (see the "Methods" section). The AB and BA MPDs have equal and opposite winding. The total winding in each moiré cell is zero, but individually the AB and BA domains have winding numbers of ± 1 2 , meaning the MPDs form a triangular network of merons and antimerons. The winding is concentrated at the domain centers and along the lines joining the AB/BA domain centers to the AA stacking configurations, and is zero along the domain walls. The magnitude of the winding is the same for strained and twisted bilayers apart from a reorientation of the axes, but the type of winding is different in each case. For strained bilayers, the merons are of Néel type, where the polarization flows into and out of the domain centers. For twisted bilayers, the merons are of Bloch type, where the polarization curls around the domain centers.

Discussion
In this work, we illustrate that, in addition to the out-of-plane polarization in layered systems like bilayer hBN, there is also an in-plane polarization as a result of a relative displacement between the layers. This phenomenon is general to all layered systems, provided that the bilayer lacks inversion symmetry. For the case of 3R-stacked hBN and similar materials (MoS 2 , etc.), there are three out-of-plane mirror planes related by C 3 rotations plus an in-plane mirror plane half-way between the layers, various combinations of which are broken as one layer slides over the other.
Our findings indicate that the polar properties of layered systems are much richer than previously thought. For untwisted bilayer hBN, the energetically favorable AB and BA stackings have zero in-plane polarization. However, knowledge of how the in-plane polarization changes during the process of vdW sliding may prove useful. For example, by measuring the change in out-of-plane polarization, through a change in out-of-plane current, it is possible to determine when vdW sliding occurs between AB and BA domains, but it is not possible to determine in which direction the sliding occurs, i.e. to which of the three neighboring domains. By measuring the change in in-plane polarization, it may be possible to distinguish between these three sliding processes, which may enhance the capacity for information processing in ferroelectric layered systems.
Ferroelectric materials have been fabricated in many different geometries, from 2D thin films and FE/PE superlattices to 1D nanowires 26,27 and nanotubes 28,29 (in fact, all nanotubes are inherently polar via flexoelectricity 30-32 ), and 0D quantum dots 33,34 . Lowerdimensional ferroelectric systems typically exhibit size-dependent transitions in which the local polarization is more complex and can exhibit vortices before the polarization eventually vanishes completely [35][36][37][38] . Soon after, it was realized these polar structures with vortices were topologically nontrivial 39 , and skyrmion-like polarization structures were then identified, for example in barium titanate (BaTiO 3 ) nanowires embedded in a matrix of strontium titanate (SrTiO 3 ) 40 . It has also been proposed that skymrions may be created by controlling domains and domain walls in ferroelectrics, where at low a Out-of-plane polarization P ⊥ (r), b change in in-plane polarization P ∥ (r) and c local topological charge q(r) for biaxially strained hBN. d out-of-plane polarization, e change in in-plane polarization, and f local topological charge for twisted hBN. Plots are shown in real space r = (r x , r y ), with with length scales of the bilayer lattice constant |a| divided by strain η and twist angle θ for strained and twisted bilayers, respectively. The black arrows represent the Moiré superlattice vectors in each case. g-i Proposal of a system in which polar topological phase transitions may be driven by an applied electric field E. A single cell of a moiré superlattice is embedded in a dielectric medium. g At zero electric field, a meron-antimeron pair forms, i.e. the topological charges in the individual domains are Q = ± 1 2 . h When a positive electric field is applied, the dielectric medium has normalized polarization P z = +1, changing the winding around the boundary of the cell. As a result, the antimeron turns into an antiskyrmion (Q → −1), and the meron vanishes (Q → 0). i When a negative field is applied, the reverse occurs: the meron turns into a skyrmion (Q → +1), and the antimeron vanishes (Q → 0). temperatures the domain walls are of Bloch type and contain an inplane polarization 41 . Ferroelectric skyrmions have recently been experimentally observed in ferroelectric/paraelectric superlattices 12,13 . In addition to skyrmions, polar merons have also been considered theoretically and signaled experimentally 14 .
A full characterization of topological polarization also provides for an interesting avenue to explore. One major conceptual problem is that, unlike other topological invariants, the description of polarization in terms of exponentially localized Wannier functions requires topologically trivial electronic bands. It may be that topological polarization is a real space analog to the topology of electronic bands in momentum space, with inversion symmetry and electric fields playing the role of time-reversal symmetry and magnetic fields. However, there is much work to be done in order to obtain a better understanding of topological polarization, and its relation to topological insulators in the presence of symmetries [42][43][44][45] as well as recently discovered multi-gap topological states due to the natural reality condition 46,47 . Nonetheless, it is evident that the topological polarization in moiré superlattices may serve as a new platform for topological physics in real materials, with great potential for the observation and control of topological phases in 2D materials.
We propose that it may be possible to drive a topological phase transition using a single moiré supercell embedded in a dielectric medium, from a single meron-antimeron pair into a skyrmion or antiskyrmion, for a fixed applied electric field (see Fig. 3g-i). Such a setup may be possible by embedding a moiré quantum dot in a dielectric material, or in an aligned bilayer in which there is a local strain or twisting around a defect, such that the bilayer is strained/ twisted inside a domain and unstrained/untwisted outside. At zero electric field, a meron-antimeron pair forms in the supercell, and the polarization in the dielectric medium is zero. When an electric field is applied, the normalized polarization in the dielectric medium is ±ẑ, which changes the winding along the boundary of the cell, adding to the winding in one domain, turning the meron/antimeron into a skyrmion/antiskyrmion, and canceling the winding in the other domain, making it topologically trivial, reminiscent of a bulk-boundary correspondence. A similar effect was also observed in periodic superlattices when the dielectric response of the system was taken into consideration. As a result, the polarization field has an in-plane component ϵ 0 ϵ ? E: everywhere in the moiré superlattice, including the domain walls. This leads to an induced winding along the domain walls, suggesting that even in periodic moiré systems, the merons/antimerons can be promoted to skyrmions/antiskyrmions, making the others topologically trivial, with an applied electric field.
It may also be possible to manipulate the meron-antimeron domain network via lattice reconstruction at different strains or twist angles. For small strains or twist angles, significant lattice relaxation can occur in order to increase the area of the more energetically stable stacking domains 24 , making the polar domains sharper 2,3 . At zero electric field, the AB and BA domains in 3R-stacked hBN relax evenly, leading to sharp triangular polar domains. However, the positions of the AA, AB, BA, and SP stackings are all preserved. As a function of strain or twist angle, this is a continuous deformation that does not break any symmetries, and the meron-antimeron network should therefore be robust against lattice relaxation at zero electric fields. When a field is applied, one type of domain will grow and the other will shrink, bending the polar domain walls. The robustness/fragility of the meron-antimeron network in response to the motion of the domain walls is not immediately clear. This goes beyond the scope of this work, however, and we leave it as a direction for future research.
In summary, we have illustrated that electronic out-of-plane and in-plane charge transfer and polarization are fundamental properties of layered systems without inversion symmetry. This will have farreaching consequences both in terms of fundamental physics in strained/twisted and commensurate layered systems, such as ferroelectricity and topology, as well as potential applications for ferroelectric-based nanodevices comprised of layers of 2D materials.

First-principles calculations
First-principles density functional theory (DFT) calculations were performed using the SIESTA 22 and ABINIT 23 codes, using PSML 48 normconserving pseudopotentials 49 , obtained from Pseudo-Dojo 50 . SIESTA employs a basis set of numerical atomic orbitals (NAOs) 22 , and double-ζ polarized (DZP) orbitals were used for all calculations. The basis sets in SIESTA were optimized by hand, following the methodology in ref. 51 . ABINIT employs a plane wave basis set, which was determined using a kinetic energy cutoff of 1000 eV. A mesh cutoff of 1200 Ry was used for the real space grid in all SIESTA calculations. A Monkhorst-Pack k-point grid 52 of 12 × 12 × 1 was used for the initial geometry relaxations, and a mesh of 18 × 18 × 1 was used to calculate the polarization. Calculations were converged until the relative changes in the Hamiltonian and density matrix were both <10 −6 . In both codes, the revPBE exchangecorrelation functional was used 53 . The C09 54,55 van der Waals correction was used in the SIESTA calculations and the vdw-DFT-D3(BJ) 56 correction was used in ABINIT. In SIESTA, when an out-of-plane electric field was applied, a dipole correction 57,58 was used in the vacuum region to prevent dipole-dipole interactions between periodic images. A dipole cutoff in slab-like systems has not been implemented in ABINIT, so although a vacuum space of 50 Å was used to separate the periodic images, the polarization is still slightly enhanced due to dipole-dipole interactions.
The top layer was translated along the unit cell diagonal over the bottom layer, which was held fixed. At each point, a geometry relaxation was performed to obtain the equilibrium layer separation, while keeping the in-plane lattice vectors fixed. The out-of-plane and in-plane polarization were then obtained by calculating the Berry phases of the Bloch states. The data were fitted to Fourier expansions which respect the C 3 rotation symmetry of bilayer hBN. It was found that both the out-of-plane and in-plane polarization were well described by the first order in the expansions, i.e. Eqs. (1) and (2). At each point along the unit cell diagonal, DFPT calculations were performed using ABINIT to calculate the dynamical charges.

Topological charge
The winding number, or topological charge, of the polarization field in configuration space is which can be mapped to real space in a strained or twisted bilayer as mentioned in the main text. Calculating the topological charge of a polarization field presents two additional complications when compared with magnetic fields. Firstly, the polarization is not of unit length and must be normalized. Secondly, there are regions in space with zero polarization (modulo a quantum of polarization), i.e. at the AA stacking configurations, near which Eq. (3) diverges. This can be avoided by calculating the topological charge following the methodology in ref. 59 . The polarization in the unit cell is discretized on a fine grid with spacing Δ. A plaquette is constructed around each grid point (see Fig. 4). The plaquettes form a grid that is offset from the original by half a grid spacing (a similar technique is used in first-principles calculations for more efficient Brillouin zone integrations 52 ). The zeros in polarization at the AA stacking configurations are thus not included in the offset grid. The local topological charge can then be defined as qðsÞ = 1 4π AðP 1 ,P 2 ,P 3 Þ + AðP 1 ,P 3 ,P 4 Þ À Á ð4Þ Article https://doi.org/10.1038/s41467-023-37337-8 where A is the signed area spanned by three points on a sphere: AðP 1 ,P 2 ,P 3 Þ = 2 arg 1 + P 1 Á P 2 + P 2 Á P 3 + P 3 Á P 1 + iP 1 Á ðP 2 × P 3 Þ À Á ð5Þ The total charge is then The total Q in the configuration space unit cell sums to zero, with the precision of around 10 −12 even for relatively coarse grids. The winding numbers of the MPDs converge to Q AB = À Q BA = 1 2 for grid spacings below Δ = 10 −4 (see Supplementary Fig. 5).
The meron-antimeron pair → (anti)skyrmion transition driven by an applied field was identified by calculating the topological charge using fixed boundary conditions outside the moiré cell rather than periodic boundary conditions. The normalized polarization in the cell is the same, and the normalized polarization outside the cell is taken to be sgn E ð Þ, where E is an applied field in the out-of-plane direction. The winding in the interior of the cell is unaffected, but an additional winding is induced along the boundary when a field is applied (see Supplementary Fig. 6). The total winding along the boundary sums to ÀsgnðEÞ. For sgnðEÞ = ± 1, the boundary contributes an additional winding of ∓ 1 2 each to the meron/antimeron, promoting one to a skymrion/antiskyrmion, and making the other topologically trivial.

Data availability
The data presented in this study were generated using free and opensource first-principles packages as described in the "Methods" section. The datasets generated during and/or analyzed during this study are available from the corresponding author upon request.

Code availability
Code used to generate the plotted polarization structures is available from the corresponding author upon request.